#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <mpi.h>
#include <unistd.h>


#define mpi_root 0

double H[3][3],b[3][3], pivote;
int i,j,k,l,cantnodos,myid,cantxnodo;
double cantidad=0;						//Cantidad de simbolos a enviar.
double cantISI=0;  						 //Cantidad de coeficientes ISI.
double SNR=0;									 //Relacion Señal-Ruido.
double *vectorISI;
double *convolucionado;
double *simbolos;
double *recibidos;

void fila_pivote(void);
void col_pivote(void);
void otros(void);
void invierte(void);
double VAGaussiana(void);

int main( int argc, char *argv[] ) {

/* Iniciamos el MPI. En cantnodos nos van a quedar el numero de nodos 
 * que metio el usuario. */	
	MPI_Init(&argc,&argv);
    MPI_Comm_size( MPI_COMM_WORLD, &cantnodos );
    MPI_Comm_rank(MPI_COMM_WORLD, &myid);
	
	double inicio=0,fin=0;
	
	
/*---Asignamos la cantidad de simbolos a enviar y la cant de coef ISI---*/
/* Solo el root lo hace. Si el usuario ingresa numeros decimales para
 * cantidad de simboos o coeficientes isi, se castea a la parte entera.*/		
	
	if(myid == mpi_root)
	{
	
	cantidad =atoi(argv[1]);
	cantISI =atoi(argv[2]);

/*------------Generacion del vector de coeficientes ISI-------------- */

	vectorISI=(double*)malloc(cantISI*sizeof(double));
	if(vectorISI==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
	/*printf("Ingrese los coeficientes ISI uno por uno: \n");
	double coef=0;
	for(i=0;i<cantISI;i++)
	{printf("%d coef: ",i+1);
	 scanf("%lf",&coef);
	 vectorISI[i]=coef;
	}*/
	
	if(cantISI==3)
	{vectorISI[0]=0.4,vectorISI[1]=0.6,vectorISI[2]=-0.2;}
	else
	{if(cantISI==5)
	 {vectorISI[0]=0.4,vectorISI[1]=0.6,vectorISI[2]=-0.2,vectorISI[3]=0.7,vectorISI[4]=-0.8;}
	 else
	 {if(cantISI==7)
	  {
	  vectorISI[0]=0.4,vectorISI[1]=0.6,vectorISI[2]=-0.2,vectorISI[3]=0.7,vectorISI[4]=-0.8,vectorISI[5]=0.1,vectorISI[6]=0.3;}
      else{exit(0);}}
	}
	 
/*--------------------------------------------------------------------*/	
	
/*-----------------------Asignamos la SNR-----------------------------*/
	SNR=atof(argv[3]);
/*--------------------------------------------------------------------*/		

	inicio= MPI_Wtime();//A partir de este momento empezamos a contar el tiempo de ejecucion.
	
	}
	MPI_Barrier(MPI_COMM_WORLD);

	
/* Cada nodo tiene un buffer de envio, y recepcion.*/	
	double *buf_recepcion,*buf_envio;								

/* --------Alocamos la memoria para el buffer de envio de cada nodo---*/	
	buf_envio=(double*)malloc(sizeof(double));
	if(buf_envio==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
	
/* Ahora el root calcula cuantos datos tiene que hacer cada nodo*/	
	if(myid == mpi_root){
	    cantxnodo=(int)cantidad/cantnodos;		
		buf_envio[0]=cantxnodo;
		buf_recepcion=(double*)malloc(cantxnodo*cantnodos*sizeof(double));
		if(buf_recepcion==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
		}
		MPI_Barrier(MPI_COMM_WORLD);
	
	
/* Mandamos a cada nodo, cuantos simbolos tiene que calcular */
	MPI_Bcast (buf_envio,1,MPI_DOUBLE,mpi_root,MPI_COMM_WORLD);
	

/* Ahora Cada nodo genera su cantidad de simbolos */
	cantxnodo=(int)buf_envio[0];//Cada nodo tiene su cantidad de simbolos.	
	free(buf_envio);
	//printf("Myid %d tengo que hacer %d ",myid,cantxnodo);
	
	simbolos=(double*)malloc(cantxnodo*cantnodos*sizeof(double));
	if(simbolos==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
	convolucionado=(double*)malloc(cantxnodo*cantnodos*sizeof(double));
	if(convolucionado==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
	recibidos=(double*)malloc(cantxnodo*cantnodos*sizeof(double));
	if(recibidos==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
	buf_envio=(double*)malloc(cantxnodo*sizeof(double));
	if(buf_envio==NULL){printf("ERROR: Out of memory\n\n");exit(0);}	

/*----------------Generacion del vector de simbolos-------------------*/
/* Antes de ingresar al ciclo for seteamos la semilla con la hora del 
 * sistema. En cada iteracion obtenemos un numero aleatorio entre 0 y 
 * 10000. Si es mayor que 5000 lo asignamos a 1 y sino a -1. De esta 
 * manera ambos simbolos son equiprobables.*/

//Cada nodo debe insertar una semilla distinta, sino se repite la secuencia.	
	srand(time(NULL)*i*getpid());													

	int simbol=0;						//Simbol tiene el num aleatorio.
	//printf("\nMyID %d\n",myid); 
	for(i=0;i<cantxnodo;i++)
	{
     simbol=rand()%10000;												
     if(simbol>5000){buf_envio[i]=1;}
	 else{buf_envio[i]=-1;}   
	
	//printf(" Sim: %4.2f ",myid,buf_envio[i]);
	}
	//printf("\nSoy %d e hice %d simbols\n",myid,i);
/*--------------------------------------------------------------------*/
	
/* Cada nodo genero su parte de los simbolos, hay que reunirlos */
	MPI_Barrier(MPI_COMM_WORLD);
	MPI_Gather(buf_envio,cantxnodo,  MPI_DOUBLE, 
				buf_recepcion,cantxnodo,  MPI_DOUBLE, 
	           	mpi_root,                  
            	MPI_COMM_WORLD);
	free(buf_envio);
/* El root reune los datos, para hacer la convolucion que se debe hacer
 *  en un procesador. En simbolos van a estar los simbolos originales a
 *  enviar*/
	
	if(myid == mpi_root)
	{ 
	//  printf("\nSoy root, hago gather\n");
	  for(i=0;i<(cantnodos*cantxnodo);i++)
	  {
	  simbolos[i]=buf_recepcion[i];
	//  printf(" %4.2f  ",simbolos[i]);
	  }
	//  printf("\nReuni %d simbolos\n",i);
	
/*-----------Convolucion del vector ISI con el de simbolos------------*/
/* En convolucionado quedan los simbolos afectados por ISI.*/
	//printf("CONVOLUCIONADO\n");
	for(i=0;i<(cantnodos*cantxnodo);i++)
	{
	 for(j=0;j<cantISI;j++)
	 {if((i-j)>=0){convolucionado[i]=convolucionado[i]+simbolos[i-j]*vectorISI[j];}
	 }
	//printf(" %4.2f ",convolucionado[i]);
	}
	//printf("\nConvolucione %d simbolos\n",i);
	free(vectorISI);
/*--------------------------------------------------------------------*/
	}
	MPI_Barrier(MPI_COMM_WORLD);


/* ----------------Paralelizacion de ruido------------------------*/
/* Paralelizamos la generacion de ruido. Primero tenemos que hacer 
 * broadcast con el SNR */
 	buf_envio=(double*)malloc(sizeof(double));
	if(buf_envio==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
	if(myid == mpi_root)
	{
	buf_envio[0]=SNR;
	}
	MPI_Barrier(MPI_COMM_WORLD);
	
	MPI_Bcast (buf_envio,1,MPI_DOUBLE,mpi_root,MPI_COMM_WORLD);
	
	SNR=buf_envio[0];		   //Todos los nodos tienen el valor de SNR.
	free(buf_envio);
/* ----Tenemos que dividir el vector de señal (convolucionado)--------*/
	if(myid == mpi_root)
	{buf_envio=(double*)malloc(cantnodos*cantxnodo*sizeof(double));
	 if(buf_envio==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
	 for(i=0;i<cantxnodo*cantnodos;i++)
		{ buf_envio[i]=convolucionado[i];}
	}
	
buf_recepcion=(double*)malloc(cantxnodo*sizeof(double));
if(buf_recepcion==NULL){printf("ERROR: Out of memory\n\n");exit(0);}

	MPI_Barrier(MPI_COMM_WORLD);
	MPI_Scatter(buf_envio,cantxnodo,  MPI_DOUBLE, 
						buf_recepcion,cantxnodo,  MPI_DOUBLE, 
	                 	mpi_root,                  
	                 	MPI_COMM_WORLD);
	
/*--------------------------------------------------------------------*/


/* ----------------Generacion y suma del ruido------------------------*/
/* Para cada simbolo vamos a calcular un nivel de ruido diferente. Cada
 * valor de ruido se calcula con VAGaussiana que entrega un numero 
 * uniformemente distribuido*/
	
	buf_envio=(double*)malloc(cantxnodo*sizeof(double));
	if(buf_envio==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
	
	double ruido=0;
	//printf("Myid: %d CON RUIDO\n",myid);
	for(i=0;i<cantxnodo;i++)
	{buf_envio[i]=buf_recepcion[i]; 
	 //printf("simbol %4.2f ",buf_envio[i]);
	 ruido = VAGaussiana();
	 buf_envio[i]=buf_envio[i]+(ruido/SNR);
	 //printf("con ruido %4.2f ",buf_envio[i]);
	}
	//printf("\n");
	free(buf_recepcion);
/*--------------------------------------------------------------------*/

/*---------Tenemos que reunir los datos con ruido---------------------*/
if(myid==mpi_root)
{buf_recepcion=(double*)malloc(cantxnodo*cantnodos*sizeof(double));
if(buf_recepcion==NULL){printf("ERROR: Out of memory\n\n");exit(0);}}
	MPI_Gather(buf_envio,cantxnodo,  MPI_DOUBLE, 
						buf_recepcion,cantxnodo,  MPI_DOUBLE, 
	                 	mpi_root,                  
	                 	MPI_COMM_WORLD);

	
/*--------------------------------------------------------------------*/
	

/* La recepcion la debe hacer el root solo (la ecualizacion). En recibidos
 * queda la señal sobre la que hay que decidir*/
 free(buf_envio);
 buf_envio=(double*)malloc(cantnodos*cantxnodo*sizeof(double));
 if(buf_envio==NULL){printf("ERROR: Out of memory\n\n");exit(0);}
 
if(myid == mpi_root)
{
		
/*----------------- Ecualizacion en la recepcion----------------------*/
/* Lo vamos a hacer mediante las matrices D=H.W
 * 0	ISI[0] ISI[1] ISI[2]	 W[0]
 * 1 =	  0	   ISI[0] ISI[1]  x	 W[1]
 * 0 	  0		 0	  ISI[0]     W[2] 
 * Entonces tenemos que encontrar W y convolucionarlo la señal que llega.
 * W= D . H^-1. Esta solucion es trivial ya que D es [0 1 0]. Se asigna
 * directamente.
 * Luego se evaluaran con el nivel de deteccion */

  H[0][0]=vectorISI[0];
  H[0][1]=vectorISI[1];
  H[0][2]=vectorISI[2];
  H[1][0]=0;
  H[1][1]=vectorISI[0];
  H[1][2]=vectorISI[1];
  H[2][0]=0;
  H[2][1]=0;
  H[2][2]=vectorISI[0];
  invierte();
/*  printf("\nMatriz inversa:\n");
 for(i=0; i<3; i++)
 {
 for(j=0; j<3; j++)
 printf("%4.2f ", H[i][j]);
 printf("\n");
 }
*/

/*---Asignamos W. Es el resultado de multiplicar H^-1 por [0 1 0]---*/
 double W[3];
 W[0]=H[1][2];
 W[1]=H[2][2];
 W[2]=H[3][2];

/*Convolucionamos. En recibidos queda la señal sobre la que hay que decidir.*/

 
 //printf("Convolucion de recepcion (con W)\n");
 for(i=0;i<(cantnodos*cantxnodo);i++)
 {convolucionado[i]=buf_recepcion[i];
  for(j=0;j<3;j++)
  {if((i-j)>=0){recibidos[i]=recibidos[i]+convolucionado[i-j]*W[j];}
  }
 // printf(" %4.2f ",recibidos[i]);
 }
 // printf("\n Recibido \n");
  for(i=0;i<cantnodos*cantxnodo;i++)
  {
  buf_envio[i]=recibidos[i];
 // printf(" %4.2f ",recibidos[i]);
  }

}//Hasta aca el root.
 MPI_Barrier(MPI_COMM_WORLD);
 
/*Hacemos broadcast de recibidos y de simbolos originales para paralelizar 
 * la deteccion y el conteo de ruido */
	
	MPI_Bcast (buf_envio,cantnodos*cantxnodo,MPI_DOUBLE,mpi_root,MPI_COMM_WORLD);
	
	for(i=0;i<cantnodos*cantxnodo;i++){	
	recibidos[i]=buf_envio[i];//Todos los nodos tienen el vector convolucionado.
	}
	//Copiamos los simbolos originales al buffer para enviar.
	if(myid==mpi_root)
	{//printf("\n Simbolos \n ");
	 for(i=0;i<cantnodos*cantxnodo;i++){	
	 buf_envio[i]=simbolos[i];
	 //printf(" %4.2f ",simbolos[i]);
	 }}
	MPI_Barrier(MPI_COMM_WORLD);
	MPI_Bcast (buf_envio,cantnodos*cantxnodo,MPI_DOUBLE,mpi_root,MPI_COMM_WORLD);
	
	for(i=0;i<cantnodos*cantxnodo;i++){	
	simbolos[i]=buf_envio[i];//Todos los nodos tienen el vector de simbolos original.
	}
/*--------------------------------------------------------------------*/

/*-------------Vamos a paralelizar la deteccion.----------------------*/ 
 int error=0;
 //printf("RECIBIDOS\n");
 for(i=0;i<cantxnodo;i++)
 {
 if(recibidos[i+(myid*cantxnodo)]>=0){recibidos[i+(myid*cantxnodo)]=1;}
  else{recibidos[i+(myid*cantxnodo)]=-1;}
//  printf(" %4.2f ",buf_envio[i]);
 //Contamos errores
 if(recibidos[i+(myid*cantxnodo)]!=simbolos[i+(myid*cantxnodo)]){error=error+1;}
 }
 
  free(buf_envio);
 buf_envio=(double*)malloc(sizeof(double));
 if(buf_envio==NULL){printf("ERROR: Out of memory\n\n");exit(0);}

 buf_envio[0]=error;
/* Por ultimo el root reune los datos del conteo de errores*/	
	
	buf_recepcion=(double*)malloc(sizeof(double));
	if(buf_recepcion==NULL){printf("ERROR: Out of memory\n\n");exit(0);}

//Hacemos reduccion
MPI_Barrier(MPI_COMM_WORLD);
	MPI_Reduce (buf_envio,buf_recepcion,1,MPI_DOUBLE,MPI_SUM,mpi_root,MPI_COMM_WORLD);

 if(myid==mpi_root){

 error=buf_recepcion[0];	 
 double BER=error*100/cantidad;
  printf("\n\n Errores: %d \n",error);
  printf("BER: %4.1f porciento\n",BER);
  fin=MPI_Wtime();
  double t_total=fin-inicio;
  printf("Tiempo de ejecucion: %f segundos\n\n",t_total);
 }
   MPI_Barrier(MPI_COMM_WORLD);
 
	free(simbolos);
	free(recibidos);
	free(convolucionado);
	free(buf_envio);
	free(buf_recepcion);
 
    MPI_Finalize();
	
return 0;
}


/*----------Funciones de Inversion de Matriz 3x3----------------------*/
void invierte()
{
for(i=0; i<3; i++)
{
j=i;
pivote=H[i][j];
b[i][j]=1/pivote;
fila_pivote();
col_pivote();
otros();
for(k=0; k<3; k++)
for(l=0; l<3; l++)
H[k][l]=b[k][l];
}	
}

void fila_pivote(void)
{
int m;
for(m=0; m<3; m++)
if(m != i)
b[i][m]=H[i][m]/pivote;
}

void col_pivote()
{
int m;
for(m=0; m<3; m++)
if(m != j)
b[m][j]=-H[m][j]/pivote;
}

void otros(void)
{
int x,y;
for(x=0 ;x<3; x++)
for(y=0; y<3; y++)
if(x!=i && y!=j)
b[x][y]=H[x][y]-(H[i][y]*H[x][j])/pivote;
}
/*--------------------------------------------------------------------*/

/*--------------Generacion de Var Aleat Gaussiana---------------------*/
double VAGaussiana(void)
{
int i, desv=1, media = 0, h = 100;
double sum=0; //Suma de h números aleatorios entre 0 y 1.

for(i=0;i<h;i++){
sum=sum+(rand()/((double)RAND_MAX));
}
return media+desv*((sum-h/2)/2.88);//2.88 = sqrt(100/12)
}
/*--------------------------------------------------------------------*/
